## In many cases, you will first need to run the code to generate the data/replicate the results before you can produce the graph

rm(list=ls())
library(ggplot2)
library(coda)
library(foreign)
library("RColorBrewer")
library(grid)



## Figure 1

statadata <- read.dta("1_elitepartychange.dta")
statadata <- statadata[which(statadata$cong>=83),]
angle <- rep(statadata$partydegree[1+2*(1:30)-2],each=2)

Measure <- factor(rep(1:2,each=60),labels=c("Avg. Angle of Cutting Line for Party Votes by Congress","Absolute Change"))
ggplotdata <- data.frame(Degrees=c(angle,c(NA,NA,statadata$partyideochange[-(1:2)])),Year=rep(statadata$year,2),Measure=Measure)

ideo1 <- ggplotdata[1:60,]
ideo2 <- ggplotdata[61:120,]


angleplot1 <- qplot(Year,Degrees,geom="step",data=ideo1) + ylab("Avg. Cutting Line Angle of Party Votes by Congress") + xlab("") + theme_bw() + theme(axis.ticks.x=element_blank(),axis.text.x=element_blank(),plot.margin=unit(c(0,.25,0,0),"cm"))
angleplot2 <- qplot(Year,Degrees,geom="step",data=ideo2) + ylab("Absolute Change") + scale_y_continuous(breaks=c(0,5,10)) + theme_bw() + theme(plot.margin=unit(c(0,.25,0,.2),"cm"))

pdf("figure1.pdf",width=8,height=6)
grid.newpage()
pushViewport(viewport(layout=grid.layout(10,1)))
vplayout <- function(x,y)
    viewport(layout.pos.row=x,layout.pos.col=y)
print(angleplot1,vp=vplayout(1:7,1))
print(angleplot2,vp=vplayout(8:10,1))
dev.off()


## Figure 2

Issue <- factor(rep(2:6,each=60),labels=c("Taxes & Management","Social Welfare","Agriculture","Civil Liberties","Foreign & Defense Policy")) 
ggplotdata <- data.frame(Year=rep(statadata$year,5),Percent=c(statadata$partyaage1*100,statadata$partyaage2*100,statadata$partyaage3*100,statadata$partyaage4*100,statadata$partyaage5*100),Issue=Issue)
mypal <- brewer.pal(8,"YlGnBu")

aageplot1 <- ggplot(ggplotdata,aes(Year,Percent)) + geom_step(aes(linetype=Issue,colour=Issue)) + scale_linetype_manual("Issue Area",values=c("Taxes & Management"=1,"Social Welfare"=2,"Agriculture"=3,"Civil Liberties"=4,"Foreign & Defense Policy"=5)) + scale_colour_manual("Issue Area",values=c("Taxes & Management"=mypal[4],"Social Welfare"=mypal[5],"Agriculture"=mypal[6],"Civil Liberties"=mypal[7],"Foreign & Defense Policy"=mypal[8])) + xlab("") + theme_bw() + theme(legend.position="top", legend.margin = unit(0, "cm"),axis.ticks.x=element_blank(),axis.text.x=element_blank(),plot.margin=unit(c(0,.25,0,0),"cm"))

ggplotdata2 <- data.frame(Year=statadata$year,Percent=statadata$partyagendachange)
aageplot2 <- qplot(Year,Percent,geom="step",data=ggplotdata2) + ylab("Overall Change") + theme_bw() + theme(plot.margin=unit(c(0,.25,0,.2),"cm"))

pdf("figure2.pdf", width=8,height=6)
grid.newpage()
pushViewport(viewport(layout=grid.layout(10,1)))
vplayout <- function(x,y)
    viewport(layout.pos.row=x,layout.pos.col=y)
print(aageplot1,vp=vplayout(1:7,1))
print(aageplot2,vp=vplayout(8:10,1))
dev.off()

## Figure 3
rm(list=ls())
statadata <- read.dta("macrocount_normally.dta")
statadata <- statadata[-(1:5),]
Macropartisanship <- c(statadata$nes_macro,statadata$gallup_macro_landline,statadata$gallup_macro_personal)*100
Quarter <- rep(seq(as.Date("1953/2/1"),as.Date("2012/11/1"),"3 months"),3)
Source <-c(rep("ANES",240),rep("Gallup Quarterly (Telephone)",240), rep("Gallup Quarterly (Personal)",240))
forggplot <- data.frame(Macropartisanship=Macropartisanship,Date=Quarter,Source=Source)

## fill in the few missing months with straight line
forggplot$Macropartisanship[623] <- (forggplot$Macropartisanship[622] + forggplot$Macropartisanship[624])/2
forggplot$Macropartisanship[627] <- (2*forggplot$Macropartisanship[626] + forggplot$Macropartisanship[629])/3
forggplot$Macropartisanship[628] <- (forggplot$Macropartisanship[626] + 2*forggplot$Macropartisanship[629])/3
forggplot$Macropartisanship[636] <- (forggplot$Macropartisanship[635] + forggplot$Macropartisanship[637])/2
forggplot$Macropartisanship[639] <- (2*forggplot$Macropartisanship[638] + forggplot$Macropartisanship[641])/3
forggplot$Macropartisanship[640] <- (forggplot$Macropartisanship[638] + 2*forggplot$Macropartisanship[641])/3
forggplot$Macropartisanship[658] <- (forggplot$Macropartisanship[657] + forggplot$Macropartisanship[659])/2


pdf("figure3.pdf",width=8,height=6)
ggplot(forggplot,aes(Quarter,Macropartisanship)) + geom_point(aes(shape=Source)) + geom_line(aes(linetype=Source)) + scale_shape_manual("Source:",values=c("ANES"=16,"Gallup Quarterly (Personal)"=26,"Gallup Quarterly (Telephone)"=26)) + scale_size_manual("Source:",values=c("ANES"=2,"Gallup Quarterly (Personal)"=1,"Gallup Quarterly (Telephone)"=1)) + scale_linetype_manual("Source:",values=c("ANES"=0,"Gallup Quarterly (Personal)"=1, "Gallup Quarterly (Telephone)"=2))   +theme_bw() + theme(legend.position=c(.8,.85),legend.title=element_text(face="bold"))
dev.off()


## need to run 2_linregresults.do & 3_MSAR_PA.R first
rm(list=ls())
load("3_MSAR_PA.Rdata")
library(coda)
library(ggplot2)


library(foreign)
mle <- read.dta("paapprovalsim.dta")
arima.mean <- apply(mle,c(2),mean)
arima.mean <- c(0,arima.mean[21:32])
arima.ci <- apply(mle,c(2),quantile,probs=c(.025,.975))
arima.ci <- cbind(c(0,0),arima.ci[,21:32])

library(grid)
class(jags.parsamples) <- "mcmc.list"

combphi1 <- c(jags.parsamples[,482][[1]],jags.parsamples[,482][[2]],jags.parsamples[,482][[3]])
combphi2 <- c(jags.parsamples[,483][[1]],jags.parsamples[,483][[2]],jags.parsamples[,483][[3]])
combphi3 <- c(jags.parsamples[,484][[1]],jags.parsamples[,484][[2]],jags.parsamples[,484][[3]])
combphi4 <- c(jags.parsamples[,485][[1]],jags.parsamples[,485][[2]],jags.parsamples[,485][[3]])


msar.omega <- (.4 + combphi2 + combphi3*.87 + combphi4*.5)/(1-combphi1)
summary(msar.omega)


approval <- c(.5,rep(.55,8),rep(.5,4))
ics <- rep(.87,13)
party <- rep(1,13)

ylag <- msar.omega
ypred <- array(NA,dim=c(12000,13))

for (i in 1:13){
   ypred[,i] <- .4 + combphi2 + combphi1*ylag + combphi3*ics[i] + combphi4*approval[i] 
   ylag <- .4 + combphi2 + combphi1*ylag + combphi3*ics[i] + combphi4*approval[i]
}

diff <- ypred-array(rep(msar.omega,13),dim=c(12000,13))
diff <- diff*100

means <- apply(diff,c(2),mean)
means[1] <- 0
ci.s <- apply(diff,c(2),quantile,probs=c(.025,.975))
ci.s
ci.s[1:2,1] <- 0

pdf("figure4.pdf",width=8,height=6)
plot(NULL,xlim=c(1,13),ylim=c(0,1.5),type="n",axes=FALSE,xlab="Quarter",ylab="Percentage Point Change in Macropartisanship")
polygon(x=c(1:13,13:1),y=c(ci.s[1,],rev(ci.s[2,])),col="gray90",lty=2,border=NA)
polygon(x=c(1:13,13:1),y=c(arima.ci[1,],rev(arima.ci[2,])),col="gray90",lty=2,border=NA)
#segments(1:13,ci.s[1,],1:13,ci.s[2,],col="gray20",lty=1)
#segments(1:13,arima.ci[1,],1:13,arima.ci[2,],col="gray20",lty=1)
abline(v=2)
abline(v=9)
lines(1:13,means,lty=1,lwd=2)
lines(1:13,arima.mean,lty=2,lwd=2)
text(4,1,expression("Approval 50% "%->%" 55%"))
text(11.2,1,expression("Approval 55% "%->%" 50%"))
legend("topright",title="Estimate",lty=c(2,1),lwd=c(2,2),c("PA","MSAR"))
axis(1)
axis(2)
dev.off()


## Figure 5

rm(list= ls())
statadata <- read.dta("4_conditionalupdatingmeasures.dta")

quarter <- seq(as.Date("1953/2/1"),as.Date("2012/11/1"),"3 months")

cumprobdis <- ifelse(statadata$avgfutureshift<.9,1,2)
cumprobdis <- ifelse(statadata$avgfutureshift<.95,cumprobdis,3)

forggplot <- data.frame(Quarter=quarter,Probability=statadata$delta,Cumprob=cumprobdis)
forggplot$Cumprob <- factor(forggplot$Cumprob,labels=c("<.9",".9-.95",">.95"))
mypal <- brewer.pal(9,"YlOrRd")
mypal <- c(mypal[3],mypal[7],mypal[9])

pdf("figure5.pdf",width=8,height=6)
ggplot(data=forggplot,aes(Quarter,Probability)) + geom_bar(stat="identity",aes(fill=Cumprob)) + scale_fill_manual("Two-Year Ahead\nProbability",values=mypal) + ylab("Probability of Mean Shift")+theme_bw() + theme(legend.position=c(.2,.85),legend.title=element_text(face="bold")) 
dev.off()

## Figure 6

my.labs <- list("Observed","Predicted",bquote(alpha[t]))
my.labs2 <- list("Observed","Predicted","Equilibrium")
mypal <- brewer.pal(3,"Set1")

 
forggplot <- data.frame(Quarter=rep(quarter,3),Macropartisanship=c(rep("Observed",240),rep("Predicted",240),rep("alpha_t",240)),Proportion=c(statadata$adjusted_gallup*100,statadata$ypred*100,statadata$ytilde*100),stringsAsFactors=FALSE)

pdf("figure6.pdf",width=8,height=6)
plot <- ggplot(forggplot,aes(Quarter,Proportion)) + geom_point(aes(shape=Macropartisanship)) + geom_line(aes(linetype=Macropartisanship,colour=Macropartisanship)) + scale_linetype_manual("Macropartisanship",breaks=c("Observed","Predicted","alpha_t"),values=c(1,0,2),labels=my.labs2) + scale_colour_manual("Macropartisanship",breaks=c("Observed","Predicted","alpha_t"),values=c(mypal[1],mypal[1],mypal[2]),labels=my.labs2) + scale_shape_manual("Macropartisanship",values=c("","o",""),breaks=c("Observed","Predicted","alpha_t"),labels=my.labs2) + ylab("Percent") + theme_bw() + theme(legend.position=c(.85,.85),legend.title=element_text(face="bold"))
plot
dev.off()

## Figure 7

forggplot <- data.frame(Quarter=rep(quarter,2),Shock=c(statadata$avgfsd,statadata$lb5fsd),Estimate=c(rep("Mean",240),rep("5th Percentile",240)))

      
mypal <- brewer.pal(9,"YlOrRd")
promplot <- ggplot(forggplot,aes(Quarter,Shock)) + geom_line(aes(linetype=Estimate,colour=Estimate)) + scale_linetype_manual("Estimate",breaks=c("Mean","5th Percentile"),values=c(2,1)) + scale_colour_manual("Estimate",breaks=c("Mean","5th Percentile"),values=c(mypal[6],mypal[9])) + ylab("Future Variance in Macropartisanship Attributed to Updating\n(Two-Year Horizon)") 
pdf("figure7.pdf",width=8,height=6)
promplot + theme_bw() + theme(legend.position=c(.8,.8),legend.title=element_text(face="bold"))
dev.off()


## Figure 8

statadata <- read.dta("avplots.dta")

year <- 1952+(1:30)*2
lb 

forggplot <- data.frame(year=rep(year,2), y = c(statadata$yagenda,statadata$yideology), x = c(statadata$agenda,statadata$ideology), variable=c(rep("Future Policy Change",30),rep("Future Alignment Change",30)), intercept=rep(0,60), slope=c(rep(.0488905,30),rep(.0212488,30)))


testplot <- ggplot(forggplot,aes(x=x,y=y,label=year))  + facet_grid(variable~.)
pdf("figure8.pdf",width=8,height=6)
testplot + geom_point() + geom_text(hjust=0,vjust=0) + geom_abline(aes(intercept=intercept,slope=slope)) + xlim(-4.1,9) + ylim(-.3,.7) + ylab("Updating\n(Deviation from Conditional Expectation)") + xlab("Deviation From Conditional Expectation") + theme_bw()
dev.off()


## Figure A.1

rm(list=ls())
statadata <- read.dta("1_gallupdataout.dta")
statadata <- statadata[-(1:129),]

Estimate <- c(statadata$gallup_macro_overall,statadata$adjusted_gallup,statadata$yhat,statadata$tbias)*100
Quarter <- rep(seq(as.Date("1984/2/1"),as.Date("2012/11/1"),"3 months"),4)
Measure <-c(rep("Observed",116),rep("Adjusted",116), rep("Smoothed Estimate",116), rep("Quarterly Estimate",116))
Facet <- c(rep("Macropartisanship",232),rep("Telephone Bias",232))

forggplot <- data.frame(Estimate,Quarter,Measure,Facet)

pdf("figureA1.pdf",width=8.5,height=7)
ggplot(forggplot,aes(Quarter,Estimate)) + geom_point(aes(shape=Measure)) + geom_line(aes(linetype=Measure)) + scale_shape_manual("Measure:",values=c("Observed"=26,"Adjusted"=26,"Quarterly Estimate"=1,"Smoothed Estimate"=26)) + scale_size_manual("Measure:",values=c("Observed"=1,"Adjusted"=1,"Quarterly Estimate"=2,"Smoothed Estimate"=1)) + scale_linetype_manual("Measure:",values=c("Observed"=2,"Adjusted"=1,"Quarterly Estimate"=0, "Smoothed Estimate"=4)) + ylab("Percentage Points") + theme_bw() + facet_grid(Facet~.,scales="free_y")  + theme(legend.title=element_text(face="bold"))
dev.off()

## Figure C.2, need to run MSAR ADL estimates first

rm(list=ls())


load("TableC3_MSAR_ADL.Rdata")
library(coda)
library(ggplot2)


library(foreign)
mle <- read.dta("adlapprovalsim.dta")
arima.mean <- apply(mle,c(2),mean)
arima.mean <- c(0,arima.mean[24:35])
arima.ci <- apply(mle,c(2),quantile,probs=c(.025,.975))
arima.ci <- cbind(c(0,0),arima.ci[,24:35])

library(grid)
class(jags.parsamples) <- "mcmc.list"
summary(jags.parsamples)


combphi1 <- c(jags.parsamples[,482][[1]],jags.parsamples[,482][[2]],jags.parsamples[,482][[3]])
combphi2 <- c(jags.parsamples[,483][[1]],jags.parsamples[,483][[2]],jags.parsamples[,483][[3]])
combphi3 <- c(jags.parsamples[,484][[1]],jags.parsamples[,484][[2]],jags.parsamples[,484][[3]])
combphi4 <- c(jags.parsamples[,485][[1]],jags.parsamples[,485][[2]],jags.parsamples[,485][[3]])
combphi5 <- c(jags.parsamples[,486][[1]],jags.parsamples[,486][[2]],jags.parsamples[,486][[3]])
combphi6 <- c(jags.parsamples[,487][[1]],jags.parsamples[,487][[2]],jags.parsamples[,487][[3]])
combphi7 <- c(jags.parsamples[,488][[1]],jags.parsamples[,488][[2]],jags.parsamples[,488][[3]])

msar.omega <- (.4 + combphi2 + combphi3*.87 + combphi4*.5 + combphi5 + combphi6*.87 + combphi7*.5)/(1-combphi1)
summary(msar.omega)


approval <- c(.5,rep(.55,8),rep(.5,4))
lagapproval <- c(.5,.5,rep(.55,8),rep(.5,3))
ics <- rep(.87,13)
party <- rep(1,13)

ylag <- msar.omega
ypred <- array(NA,dim=c(12000,13))

for (i in 1:13){
   ypred[,i] <- .4 + combphi2 + combphi1*ylag + combphi3*ics[i] + combphi4*approval[i] + combphi5 + combphi6*ics[i] + combphi7*lagapproval[i]
   ylag <- .4 + combphi2 + combphi1*ylag + combphi3*ics[i] + combphi4*approval[i] + combphi5 + combphi6*ics[i] + combphi7*lagapproval[i]
}

diff <- ypred-array(rep(msar.omega,13),dim=c(12000,13))
diff <- diff*100

means <- apply(diff,c(2),mean)
means[1] <- 0
ci.s <- apply(diff,c(2),quantile,probs=c(.025,.975))
ci.s
ci.s[1:2,1] <- 0

pdf("figureC2.pdf",width=8,height=6)
plot(NULL,xlim=c(1,13),ylim=c(0,1.5),type="n",axes=FALSE,xlab="Quarter",ylab="Percentage Point Change in Macropartisanship")
abline(v=2)
abline(v=9)
lines(1:13,means,lty=1,lwd=2)
lines(1:13,arima.mean,lty=2,lwd=2)
text(4,1,expression("Approval 50% "%->%" 55%"))
text(11.2,1,expression("Approval 55% "%->%" 50%"))
legend("topright",title="Estimate",lty=c(2,1),lwd=c(2,2),c("ADL","MSAR-ADL"))
axis(1)
axis(2)
dev.off()


## Codes and data to produce Figures C3 and C4 are in TableC4_statespace_resultsandfigures.R
